#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern double cvt_density();
extern void mincalc_flux();

extern double curvature_indicator[3];

/* second order center difference - extrapolate at the points closest
   but not on the boundary */

#define d_buf(ix,iy) gdt(dnum,(ix),(iy))

/* compute the curvature using central differencing.
   input ix,iy,dnum - level set identifier to compute curvature of.
   output - distance_dd (value of sqrt(sqr(dx)+sqr(dy)))
            dx,dy,out_distance - values at ix,iy
            return value of curvature.
*/
double compute_curvature(ix,iy,distance_dd,dx,dy,out_distance,dnum)
int ix,iy;
int dnum;
double *distance_dd,*dx,*dy,*out_distance;
{
double rplus,rminus,rmain;
double dx1[2][2],dy1[2][2],normd[2][2];
double dxx,dyy;
double curvature;
int m,n;
logical zero_norm;
double curv1,curv2,distval;

  *out_distance=d_buf(ix,iy);

  if (use_axisymmetric)
   {
    rmain=x[ix];
    rplus=rmain+h/2.0;
    rminus=rplus-h;
   }
  else
   {
    rmain=rplus=rminus=1.0;
   }

#define dval(incx,incy) (d_buf(ix+(incx),iy+(incy)) )

  dx1[0][0]=(dval(0,0)+dval(0,-1)-dval(-1,0)-dval(-1,-1))/(2.0*h);
  dx1[0][1]=(dval(0,1)+dval(0,0)-dval(-1,1)-dval(-1,0))/(2.0*h);
  dx1[1][0]=(dval(1,0)+dval(1,-1)-dval(0,0)-dval(0,-1))/(2.0*h);
  dx1[1][1]=(dval(1,1)+dval(1,0)-dval(0,1)-dval(0,0))/(2.0*h);

  dy1[0][0]=(dval(0,0)-dval(0,-1)+dval(-1,0)-dval(-1,-1))/(2.0*h);
  dy1[0][1]=(dval(0,1)-dval(0,0)+dval(-1,1)-dval(-1,0))/(2.0*h);
  dy1[1][0]=(dval(1,0)-dval(1,-1)+dval(0,0)-dval(0,-1))/(2.0*h);
  dy1[1][1]=(dval(1,1)-dval(1,0)+dval(0,1)-dval(0,0))/(2.0*h);

  zero_norm=FALSE;
  for (m=0;m<2;m++)
   for (n=0;n<2;n++)
    {
     normd[m][n]=sqrt(sqr(dx1[m][n])+sqr(dy1[m][n]));
     if (normd[m][n]<=almost_same)
      zero_norm=TRUE;
    }

  *dx=(dx1[0][0]+dx1[0][1]+dx1[1][1]+dx1[1][0])/4.0;
  *dy=(dy1[0][0]+dy1[0][1]+dy1[1][1]+dy1[1][0])/4.0;
 
  *distance_dd=sqrt(sqr(*dx)+sqr(*dy));
  if (*distance_dd<=almost_same)
   zero_norm=TRUE;

  if (! zero_norm)
   {
    dxx=(rplus*(dx1[1][1]/normd[1][1]+dx1[1][0]/normd[1][0])-
         rminus*(dx1[0][1]/normd[0][1]+dx1[0][0]/normd[0][0]))/(2*rmain*h);
    dyy=((dy1[1][1]/normd[1][1]+dy1[0][1]/normd[0][1])-
         (dy1[1][0]/normd[1][0]+dy1[0][0]/normd[0][0]))/(2*h);

    curvature=(dxx+dyy);   /* curvature is (1/R1+1/R2) */
    /* in 2d R2=infinity so the curvature is just 1/R1 */

    if (curvature_test)
     {
      if (use_axisymmetric)
       {
        dxx=((dx1[1][1]/normd[1][1]+dx1[1][0]/normd[1][0])-
             (dx1[0][1]/normd[0][1]+dx1[0][0]/normd[0][0]))/(2*h);
        curv1=dxx+dyy;
/* not good  curv2=((dx1[1][1]/normd[1][1]+dx1[1][0]/normd[1][0])+
               (dx1[0][1]/normd[0][1]+dx1[0][0]/normd[0][0]))/(4.0*rmain);
*/
        curv2=curvature-curv1;
       }
      else
       {
        curv1=curvature;
        curv2=0.0;
       }
      distval=dval(0,0);
      if (1.0-distval*curv1<=0.0) curv1=2.0/h;
      else curv1=curv1/(1.0-distval*curv1);
      if (1.0-distval*curv2<=0.0) curv2=2.0/h;
      else curv2=curv2/(1.0-distval*curv2);
      if (curv1>2.0/h) curv1=2.0/h;
      if (curv2>2.0/h) curv2=2.0/h;
      curvature=curv1+curv2;
     }   /* curvature_test */


    return curvature;
   }   /* above if !zero_norm */
  else
   return 0.0;
}   /* compute_curvature */

/* compute H(phi) where delta = d(H)/d(phi)
   input in_spread,in_distance
*/
double compute_heaviside_general(in_distance,in_spread)
double in_distance,in_spread;
{
double cutoff;

  cutoff = in_spread*h;

  if (in_distance>=cutoff)
    return 0.5;
  else if (in_distance<=-cutoff)
    return -0.5;
  else
    return in_distance/(2.0*cutoff)+sin(pi*in_distance/cutoff)/(2.0*pi);
}   /* compute_heaviside_general */

double compute_heaviside(in_distance)
double in_distance;
{
 return compute_heaviside_general(in_distance,density_spread);
}

double compute_delta_general(in_distance,in_spread)
double in_distance,in_spread;
{
double cutoff;

  cutoff = in_spread*h;
  if (fabs(in_distance)>=cutoff)
    return 0.0;

  /*return (1.0-fabs(in_distance)/cutoff)/cutoff;*/
  return (1.0+cos(pi*in_distance/cutoff))/(2.0*cutoff);
}  /* compute_delta_general */


/* compute delta(phi)
   input density_spread,in_distance
*/
double compute_delta(in_distance)
double in_distance;
{
  return compute_delta_general(in_distance,density_spread);
}  /* compute_delta */

double compute_delta_prime(in_distance,in_spread)
double in_distance,in_spread;
{
double cutoff;

  cutoff = in_spread*h;
  if (fabs(in_distance)>=cutoff)
    return 0.0;

  return -(pi/(2.0*sqr(cutoff)))*sin(pi*in_distance/cutoff);
}


double compute_heaviside_nozero(in_distance,in_spread)
double in_distance,in_spread;
{
double cutoff,temp;

  temp=sqr(1-(1.0*h/(xmax-xmin)));
  cutoff=sqr(in_spread*h)*(1.0-temp)/temp;
  return 0.5*(1.0+in_distance/(sqrt(sqr(in_distance)+cutoff)) );
}   /* compute_heaviside_general */


double compute_delta_nozero(in_distance,in_spread)
double in_distance,in_spread;
{
double cutoff,temp;

  temp=sqr(1-(1.0*h/(xmax-xmin)));
  cutoff=sqr(in_spread*h)*(1.0-temp)/temp;

  temp=sqrt(sqr(in_distance)+cutoff);
  return cutoff/(2.0*sqr(temp)*temp);
}  /* compute_delta_general */
